#include "fft_radix13.h"

#include "fft_common_use.h"

vect_complex_d fft_radix_13(vect_complex_d& f)
{
    int N = f.size();
    vect_complex_d F(N);
    if (N == 13)
    {
        F[0] = f[0] +  f[1] + f[2] + f[3] + f[4] + f[5] + f[6] + f[7] + f[8] + f[9] + f[10] + f[11] + f[12];
        F[1] = f[0] + complex_d(0.8854560256532099, -0.4647231720437685) * f[1] + complex_d(0.5680647467311559, -0.8229838658936564) * f[2] + complex_d(0.120536680255323, -0.992708874098054) * f[3] + complex_d(-0.35460488704253545, -0.9350162426854148) * f[4] + complex_d(-0.7485107481711012, -0.6631226582407952) * f[5] + complex_d(-0.970941817426052, -0.23931566428755768) * f[6] + complex_d(-0.9709418174260521, 0.23931566428755743) * f[7] + complex_d(-0.7485107481711013, 0.663122658240795) * f[8] + complex_d(-0.3546048870425359, 0.9350162426854147) * f[9] + complex_d(0.1205366802553232, 0.992708874098054) * f[10] + complex_d(0.5680647467311548, 0.822983865893657) * f[11] + complex_d(0.88545602565321, 0.4647231720437684) * f[12];
        F[2] = f[0] + complex_d(0.5680647467311559, -0.8229838658936564) * f[1] + complex_d(-0.35460488704253545, -0.9350162426854148) * f[2] + complex_d(-0.970941817426052, -0.23931566428755768) * f[3] + complex_d(-0.7485107481711013, 0.663122658240795) * f[4] + complex_d(0.1205366802553232, 0.992708874098054) * f[5] + complex_d(0.88545602565321, 0.4647231720437684) * f[6] + complex_d(0.8854560256532099, -0.4647231720437685) * f[7] + complex_d(0.120536680255323, -0.992708874098054) * f[8] + complex_d(-0.7485107481711012, -0.6631226582407952) * f[9] + complex_d(-0.9709418174260521, 0.23931566428755743) * f[10] + complex_d(-0.3546048870425359, 0.9350162426854147) * f[11] + complex_d(0.5680647467311548, 0.822983865893657) * f[12];
        F[3] = f[0] + complex_d(0.120536680255323, -0.992708874098054) * f[1] + complex_d(-0.970941817426052, -0.23931566428755768) * f[2] + complex_d(-0.3546048870425359, 0.9350162426854147) * f[3] + complex_d(0.88545602565321, 0.4647231720437684) * f[4] + complex_d(0.5680647467311559, -0.8229838658936564) * f[5] + complex_d(-0.7485107481711012, -0.6631226582407952) * f[6] + complex_d(-0.7485107481711013, 0.663122658240795) * f[7] + complex_d(0.5680647467311548, 0.822983865893657) * f[8] + complex_d(0.8854560256532099, -0.4647231720437685) * f[9] + complex_d(-0.35460488704253545, -0.9350162426854148) * f[10] + complex_d(-0.9709418174260521, 0.23931566428755743) * f[11] + complex_d(0.1205366802553232, 0.992708874098054) * f[12];
        F[4] = f[0] + complex_d(-0.35460488704253545, -0.9350162426854148) * f[1] + complex_d(-0.7485107481711013, 0.663122658240795) * f[2] + complex_d(0.88545602565321, 0.4647231720437684) * f[3] + complex_d(0.120536680255323, -0.992708874098054) * f[4] + complex_d(-0.9709418174260521, 0.23931566428755743) * f[5] + complex_d(0.5680647467311548, 0.822983865893657) * f[6] + complex_d(0.5680647467311559, -0.8229838658936564) * f[7] + complex_d(-0.970941817426052, -0.23931566428755768) * f[8] + complex_d(0.1205366802553232, 0.992708874098054) * f[9] + complex_d(0.8854560256532099, -0.4647231720437685) * f[10] + complex_d(-0.7485107481711012, -0.6631226582407952) * f[11] + complex_d(-0.3546048870425359, 0.9350162426854147) * f[12];
        F[5] = f[0] + complex_d(-0.7485107481711012, -0.6631226582407952) * f[1] + complex_d(0.1205366802553232, 0.992708874098054) * f[2] + complex_d(0.5680647467311559, -0.8229838658936564) * f[3] + complex_d(-0.9709418174260521, 0.23931566428755743) * f[4] + complex_d(0.88545602565321, 0.4647231720437684) * f[5] + complex_d(-0.35460488704253545, -0.9350162426854148) * f[6] + complex_d(-0.3546048870425359, 0.9350162426854147) * f[7] + complex_d(0.8854560256532099, -0.4647231720437685) * f[8] + complex_d(-0.970941817426052, -0.23931566428755768) * f[9] + complex_d(0.5680647467311548, 0.822983865893657) * f[10] + complex_d(0.120536680255323, -0.992708874098054) * f[11] + complex_d(-0.7485107481711013, 0.663122658240795) * f[12];
        F[6] = f[0] + complex_d(-0.970941817426052, -0.23931566428755768) * f[1] + complex_d(0.88545602565321, 0.4647231720437684) * f[2] + complex_d(-0.7485107481711012, -0.6631226582407952) * f[3] + complex_d(0.5680647467311548, 0.822983865893657) * f[4] + complex_d(-0.35460488704253545, -0.9350162426854148) * f[5] + complex_d(0.1205366802553232, 0.992708874098054) * f[6] + complex_d(0.120536680255323, -0.992708874098054) * f[7] + complex_d(-0.3546048870425359, 0.9350162426854147) * f[8] + complex_d(0.5680647467311559, -0.8229838658936564) * f[9] + complex_d(-0.7485107481711013, 0.663122658240795) * f[10] + complex_d(0.8854560256532099, -0.4647231720437685) * f[11] + complex_d(-0.9709418174260521, 0.23931566428755743) * f[12];
        F[7] = f[0] + complex_d(-0.9709418174260521, 0.23931566428755743) * f[1] + complex_d(0.8854560256532099, -0.4647231720437685) * f[2] + complex_d(-0.7485107481711013, 0.663122658240795) * f[3] + complex_d(0.5680647467311559, -0.8229838658936564) * f[4] + complex_d(-0.3546048870425359, 0.9350162426854147) * f[5] + complex_d(0.120536680255323, -0.992708874098054) * f[6] + complex_d(0.1205366802553232, 0.992708874098054) * f[7] + complex_d(-0.35460488704253545, -0.9350162426854148) * f[8] + complex_d(0.5680647467311548, 0.822983865893657) * f[9] + complex_d(-0.7485107481711012, -0.6631226582407952) * f[10] + complex_d(0.88545602565321, 0.4647231720437684) * f[11] + complex_d(-0.970941817426052, -0.23931566428755768) * f[12];
        F[8] = f[0] + complex_d(-0.7485107481711013, 0.663122658240795) * f[1] + complex_d(0.120536680255323, -0.992708874098054) * f[2] + complex_d(0.5680647467311548, 0.822983865893657) * f[3] + complex_d(-0.970941817426052, -0.23931566428755768) * f[4] + complex_d(0.8854560256532099, -0.4647231720437685) * f[5] + complex_d(-0.3546048870425359, 0.9350162426854147) * f[6] + complex_d(-0.35460488704253545, -0.9350162426854148) * f[7] + complex_d(0.88545602565321, 0.4647231720437684) * f[8] + complex_d(-0.9709418174260521, 0.23931566428755743) * f[9] + complex_d(0.5680647467311559, -0.8229838658936564) * f[10] + complex_d(0.1205366802553232, 0.992708874098054) * f[11] + complex_d(-0.7485107481711012, -0.6631226582407952) * f[12];
        F[9] = f[0] + complex_d(-0.3546048870425359, 0.9350162426854147) * f[1] + complex_d(-0.7485107481711012, -0.6631226582407952) * f[2] + complex_d(0.8854560256532099, -0.4647231720437685) * f[3] + complex_d(0.1205366802553232, 0.992708874098054) * f[4] + complex_d(-0.970941817426052, -0.23931566428755768) * f[5] + complex_d(0.5680647467311559, -0.8229838658936564) * f[6] + complex_d(0.5680647467311548, 0.822983865893657) * f[7] + complex_d(-0.9709418174260521, 0.23931566428755743) * f[8] + complex_d(0.120536680255323, -0.992708874098054) * f[9] + complex_d(0.88545602565321, 0.4647231720437684) * f[10] + complex_d(-0.7485107481711013, 0.663122658240795) * f[11] + complex_d(-0.35460488704253545, -0.9350162426854148) * f[12];
        F[10] = f[0] + complex_d(0.1205366802553232, 0.992708874098054) * f[1] + complex_d(-0.9709418174260521, 0.23931566428755743) * f[2] + complex_d(-0.35460488704253545, -0.9350162426854148) * f[3] + complex_d(0.8854560256532099, -0.4647231720437685) * f[4] + complex_d(0.5680647467311548, 0.822983865893657) * f[5] + complex_d(-0.7485107481711013, 0.663122658240795) * f[6] + complex_d(-0.7485107481711012, -0.6631226582407952) * f[7] + complex_d(0.5680647467311559, -0.8229838658936564) * f[8] + complex_d(0.88545602565321, 0.4647231720437684) * f[9] + complex_d(-0.3546048870425359, 0.9350162426854147) * f[10] + complex_d(-0.970941817426052, -0.23931566428755768) * f[11] + complex_d(0.120536680255323, -0.992708874098054) * f[12];
        F[11] = f[0] + complex_d(0.5680647467311548, 0.822983865893657) * f[1] + complex_d(-0.3546048870425359, 0.9350162426854147) * f[2] + complex_d(-0.9709418174260521, 0.23931566428755743) * f[3] + complex_d(-0.7485107481711012, -0.6631226582407952) * f[4] + complex_d(0.120536680255323, -0.992708874098054) * f[5] + complex_d(0.8854560256532099, -0.4647231720437685) * f[6] + complex_d(0.88545602565321, 0.4647231720437684) * f[7] + complex_d(0.1205366802553232, 0.992708874098054) * f[8] + complex_d(-0.7485107481711013, 0.663122658240795) * f[9] + complex_d(-0.970941817426052, -0.23931566428755768) * f[10] + complex_d(-0.35460488704253545, -0.9350162426854148) * f[11] + complex_d(0.5680647467311559, -0.8229838658936564) * f[12];
        F[12] = f[0] + complex_d(0.88545602565321, 0.4647231720437684) * f[1] + complex_d(0.5680647467311548, 0.822983865893657) * f[2] + complex_d(0.1205366802553232, 0.992708874098054) * f[3] + complex_d(-0.3546048870425359, 0.9350162426854147) * f[4] + complex_d(-0.7485107481711013, 0.663122658240795) * f[5] + complex_d(-0.9709418174260521, 0.23931566428755743) * f[6] + complex_d(-0.970941817426052, -0.23931566428755768) * f[7] + complex_d(-0.7485107481711012, -0.6631226582407952) * f[8] + complex_d(-0.35460488704253545, -0.9350162426854148) * f[9] + complex_d(0.120536680255323, -0.992708874098054) * f[10] + complex_d(0.5680647467311559, -0.8229838658936564) * f[11] + complex_d(0.8854560256532099, -0.4647231720437685) * f[12];
        return F;
    }
    int limit = N / 13;
    vect_complex_d f_1(limit), f_2(limit), f_3(limit),
        f_4(limit), f_5(limit), f_6(limit), f_7(limit),
        f_8(limit), f_9(limit), f_10(limit), f_11(limit),
        f_12(limit), f_13(limit);
    for (int r = 0; r < limit; r++)
    {
        f_1[r] = f[13 * r];
        f_2[r] = f[13 * r + 1];
        f_3[r] = f[13 * r + 2];
        f_4[r] = f[13 * r + 3];
        f_5[r] = f[13 * r + 4];
        f_6[r] = f[13 * r + 5];
        f_7[r] = f[13 * r + 6];
        f_8[r] = f[13 * r + 7];
        f_9[r] = f[13 * r + 8];
        f_10[r] = f[13 * r + 9];
        f_11[r] = f[13 * r + 10];
        f_12[r] = f[13 * r + 11];
        f_13[r] = f[13 * r + 12];
    }
    vect_complex_d F_1, F_2, F_3,
        F_4, F_5, F_6, F_7,
        F_8, F_9, F_10, F_11,
        F_12, F_13;
    F_1 = fft_radix_13(f_1);
    F_2 = fft_radix_13(f_2);
    F_3 = fft_radix_13(f_3);
    F_4 = fft_radix_13(f_4);
    F_5 = fft_radix_13(f_5);
    F_6 = fft_radix_13(f_6);
    F_7 = fft_radix_13(f_7);
    F_8 = fft_radix_13(f_8);
    F_9 = fft_radix_13(f_9);
    F_10 = fft_radix_13(f_10);
    F_11 = fft_radix_13(f_11);
    F_12 = fft_radix_13(f_12);
    F_13 = fft_radix_13(f_13);
    for (int k = 0; k < limit; k++)
    {
        complex_d W1(calculate_W(N, k));
        complex_d W2(calculate_W(N, 2 * k));
        complex_d W3(calculate_W(N, 3 * k));
        complex_d W4(calculate_W(N, 4 * k));
        complex_d W5(calculate_W(N, 5 * k));
        complex_d W6(calculate_W(N, 6 * k));
        complex_d W7(calculate_W(N, 7 * k));
        complex_d W8(calculate_W(N, 8 * k));
        complex_d W9(calculate_W(N, 9 * k));
        complex_d W10(calculate_W(N, 10 * k));
        complex_d W11(calculate_W(N, 11 * k));
        complex_d W12(calculate_W(N, 12 * k));
        F[k] = F_1[k] + W1 * F_2[k] + W2 * F_3[k] + W3 * F_4[k] + W4 * F_5[k] + W5 * F_6[k] + W6 * F_7[k] + W7 * F_8[k] + W8 * F_9[k] + W9 * F_10[k] + W10 * F_11[k] + W11 * F_12[k] + W12 * F_13[k];
        F[k + limit] = F_1[k] + w13_1 * W1 * F_2[k] + w13_2 * W2 * F_3[k] + w13_3 * W3 * F_4[k] + w13_4 * W4 * F_5[k] + w13_5 * W5 * F_6[k] + w13_6 * W6 * F_7[k] + w13_7 * W7 * F_8[k] + w13_8 * W8 * F_9[k] + w13_9 * W9 * F_10[k] + w13_10 * W10 * F_11[k] + w13_11 * W11 * F_12[k] + w13_12 * W12 * F_13[k];
        F[k + 2 * limit] = F_1[k] + w13_2 * W1 * F_2[k] + w13_4 * W2 * F_3[k] + w13_6 * W3 * F_4[k] + w13_8 * W4 * F_5[k] + w13_10 * W5 * F_6[k] + w13_12 * W6 * F_7[k] + w13_1 * W7 * F_8[k] + w13_3 * W8 * F_9[k] + w13_5 * W9 * F_10[k] + w13_7 * W10 * F_11[k] + w13_9 * W11 * F_12[k] + w13_11 * W12 * F_13[k];
        F[k + 3 * limit] = F_1[k] + w13_3 * W1 * F_2[k] + w13_6 * W2 * F_3[k] + w13_9 * W3 * F_4[k] + w13_12 * W4 * F_5[k] + w13_2 * W5 * F_6[k] + w13_5 * W6 * F_7[k] + w13_8 * W7 * F_8[k] + w13_11 * W8 * F_9[k] + w13_1 * W9 * F_10[k] + w13_4 * W10 * F_11[k] + w13_7 * W11 * F_12[k] + w13_10 * W12 * F_13[k];
        F[k + 4 * limit] = F_1[k] + w13_4 * W1 * F_2[k] + w13_8 * W2 * F_3[k] + w13_12 * W3 * F_4[k] + w13_3 * W4 * F_5[k] + w13_7 * W5 * F_6[k] + w13_11 * W6 * F_7[k] + w13_2 * W7 * F_8[k] + w13_6 * W8 * F_9[k] + w13_10 * W9 * F_10[k] + w13_1 * W10 * F_11[k] + w13_5 * W11 * F_12[k] + w13_9 * W12 * F_13[k];
        F[k + 5 * limit] = F_1[k] + w13_5 * W1 * F_2[k] + w13_10 * W2 * F_3[k] + w13_2 * W3 * F_4[k] + w13_7 * W4 * F_5[k] + w13_12 * W5 * F_6[k] + w13_4 * W6 * F_7[k] + w13_9 * W7 * F_8[k] + w13_1 * W8 * F_9[k] + w13_6 * W9 * F_10[k] + w13_11 * W10 * F_11[k] + w13_3 * W11 * F_12[k] + w13_8 * W12 * F_13[k];
        F[k + 6 * limit] = F_1[k] + w13_6 * W1 * F_2[k] + w13_12 * W2 * F_3[k] + w13_5 * W3 * F_4[k] + w13_11 * W4 * F_5[k] + w13_4 * W5 * F_6[k] + w13_10 * W6 * F_7[k] + w13_3 * W7 * F_8[k] + w13_9 * W8 * F_9[k] + w13_2 * W9 * F_10[k] + w13_8 * W10 * F_11[k] + w13_1 * W11 * F_12[k] + w13_7 * W12 * F_13[k];
        F[k + 7 * limit] = F_1[k] + w13_7 * W1 * F_2[k] + w13_1 * W2 * F_3[k] + w13_8 * W3 * F_4[k] + w13_2 * W4 * F_5[k] + w13_9 * W5 * F_6[k] + w13_3 * W6 * F_7[k] + w13_10 * W7 * F_8[k] + w13_4 * W8 * F_9[k] + w13_11 * W9 * F_10[k] + w13_5 * W10 * F_11[k] + w13_12 * W11 * F_12[k] + w13_6 * W12 * F_13[k];
        F[k + 8 * limit] = F_1[k] + w13_8 * W1 * F_2[k] + w13_3 * W2 * F_3[k] + w13_11 * W3 * F_4[k] + w13_6 * W4 * F_5[k] + w13_1 * W5 * F_6[k] + w13_9 * W6 * F_7[k] + w13_4 * W7 * F_8[k] + w13_12 * W8 * F_9[k] + w13_7 * W9 * F_10[k] + w13_2 * W10 * F_11[k] + w13_10 * W11 * F_12[k] + w13_5 * W12 * F_13[k];
        F[k + 9 * limit] = F_1[k] + w13_9 * W1 * F_2[k] + w13_5 * W2 * F_3[k] + w13_1 * W3 * F_4[k] + w13_10 * W4 * F_5[k] + w13_6 * W5 * F_6[k] + w13_2 * W6 * F_7[k] + w13_11 * W7 * F_8[k] + w13_7 * W8 * F_9[k] + w13_3 * W9 * F_10[k] + w13_12 * W10 * F_11[k] + w13_8 * W11 * F_12[k] + w13_4 * W12 * F_13[k];
        F[k + 10 * limit] = F_1[k] + w13_10 * W1 * F_2[k] + w13_7 * W2 * F_3[k] + w13_4 * W3 * F_4[k] + w13_1 * W4 * F_5[k] + w13_11 * W5 * F_6[k] + w13_8 * W6 * F_7[k] + w13_5 * W7 * F_8[k] + w13_2 * W8 * F_9[k] + w13_12 * W9 * F_10[k] + w13_9 * W10 * F_11[k] + w13_6 * W11 * F_12[k] + w13_3 * W12 * F_13[k];
        F[k + 11 * limit] = F_1[k] + w13_11 * W1 * F_2[k] + w13_9 * W2 * F_3[k] + w13_7 * W3 * F_4[k] + w13_5 * W4 * F_5[k] + w13_3 * W5 * F_6[k] + w13_1 * W6 * F_7[k] + w13_12 * W7 * F_8[k] + w13_10 * W8 * F_9[k] + w13_8 * W9 * F_10[k] + w13_6 * W10 * F_11[k] + w13_4 * W11 * F_12[k] + w13_2 * W12 * F_13[k];
        F[k + 12 * limit] = F_1[k] + w13_12 * W1 * F_2[k] + w13_11 * W2 * F_3[k] + w13_10 * W3 * F_4[k] + w13_9 * W4 * F_5[k] + w13_8 * W5 * F_6[k] + w13_7 * W6 * F_7[k] + w13_6 * W7 * F_8[k] + w13_5 * W8 * F_9[k] + w13_4 * W9 * F_10[k] + w13_3 * W10 * F_11[k] + w13_2 * W11 * F_12[k] + w13_1 * W12 * F_13[k];
    }
    return F;
}
